/**                   GNU GENERAL PUBLIC LICENSE
 *                       Version 3, 29 June 2007
 *
 * Copyright (C) 2007 Free Software Foundation, Inc. <http://fsf.org/>
 * Everyone is permitted to copy and distribute verbatim copies
 * of this license document, but changing it is not allowed.
 * See License.txt file that comes with this distribution
 *
 * Description: 
 * Equality solver. Implements Gaussian elimination.
 * Return system of equalities without the eliminated variable.
 * The size of the output matrix is equal to the size of the resulting matrix
 *
 * Author: Vitaly German <vitaly.german@gmail.com>, (C) 2006-2007
 * 
 **/

/** Elimination of all inequlities by substitution
 *
 * Solves system of inequalities by substitution of bounds
 * Input system assumed to be in form Ax >= b (InqSolverSubst)
 * or Ax <= b (InqSolverSubstLEq)
 *  
 **/

#pragma once
#include <iostream>
#include <set>
#include <stack>

#include <algorithm>

#include <eliminate_eq.hpp>
#include <boost/shared_ptr.hpp>

template <class MX, class V>
class InqSolverSubst;


struct InqSolverResult{
    InqSolverResult():consistent(true),  nbrunches(0), nsolved(0){}
    InqSolverResult(bool cons, long long nbr, long long nsol):consistent(cons),  nbrunches(nbr), nsolved(nsol)
    {}
    
    bool consistent;
    long long nbrunches;
    long long nsolved;
};


typedef  std::set<size_t> Vars;

/// I/O 
template <class T>
std::ostream & operator<< (std::ostream & o, const std::set<T> &v){
    o<<"["<<v.size()<<"]: set (";
    std::copy(v.begin(), v.end(), std::ostream_iterator<T>(o, " "));
    o<<")";
    return o;
}

template <class MX, class V>
class Env {
    MX inq_;
    V cnt_;
    Vars var_;
    size_t line_;
public:
    typedef boost::shared_ptr< Env<MX, V> > ptr;
    Env(MX &inq,  V & cnt):inq_(inq), cnt_(cnt), line_(0){}
    Env(MX &inq,  V & cnt, Vars & var, size_t line):inq_(inq), cnt_(cnt), var_(var), line_(line) {}
    Env(const Env & o):inq_ (o.inq_), cnt_(o.cnt_), var_ (o.var_), line_(o.line_){}
    const Env & operator=(const Env &o){
    	if(this != &o){
    	    inq_ = o.inq_;
            cnt_ = o.cnt_;
            var_ = o.var_;
            line_ = o.line_;
    	}
    	return *this;
    }
    MX& inq() {return inq_;};
    V& cnt() {return cnt_;}
    Vars& var() {return var_;}
    size_t& line() {return line_;}
};

template <class MX>
size_t nlive(const MX &o){
    typename MX::value_type zero(0);
    size_t clive = 0;
    for(size_t i=0; i<o.size1();++i)
    	for(size_t j=0; j<o.size2();++j)
    	    if (o(i, j)!= zero) 
    	    	clive++;
    return clive;
}

/// Solves Ax >=b problem
template <class MX, class V>
class InqSolverSubst {

    typename Env< MX, V>::ptr env;
    std::stack< typename Env< MX, V >::ptr > env_stack_;

public:

    /// 1:-1 unknown, 0- no solution 1- has solution
    /// 2: complexity (shall be very long int here)

    InqSolverSubst(MX & inq, V & c):
        env(typename Env< MX, V >::ptr(new Env< MX, V >(inq, c))) {}
    InqSolverResult Solve();
    virtual ~InqSolverSubst() {}

private:

    bool can_solve();
    Vars build_vars() const;
    /// find next inequality to solve
    /// if returns false, then there are no equalities to solve
    /// remove variable (row) from vars if all zero found in the given line
    bool find_inequ(Vars & vars, size_t & line, size_t & row) const;

    /// returns true if the is one more solution
    bool has_next(size_t & line, const size_t & row) const;
protected:
    virtual bool check_solution() const;
    const V& cnt() const {
        return env->cnt();
    }
    const MX& inq() const {
        return env->inq();
    }
};

/// Solves Ax <=b problem
template <class MX, class V >
class InqSolverSubstLEq: public InqSolverSubst<MX, V > {
public:
    InqSolverSubstLEq(MX & inq, V & c):InqSolverSubst<MX, V >(inq, c) {}
protected:
    bool check_solution() const;
};

///=================== IMPLEMENTATION ====================
template <class MX, class V >
InqSolverResult InqSolverSubst<MX, V >::Solve() {

    /// Algorithm:
    /// 1. Build set of non empty columns (variables)
    /// 2. Calculate position to eliminate (line, var)
    /// 3. if other solutions are possible, the line number is saved to new_line, a new
    /// environment (new_env) is created from the current one and 
    ///    and line varible is updated to the found one (new_env.line = new_line)
    /// 4. push the new_env to the stack
    /// 5. substitute the variable at position (line, var)
    /// 6. If all variables was elimianted check if vector of constants has
    ///    values more then 0. In this case system does not have solutions
    /// 7. If solution found - return.
    /// 8. if stack is not empty get the environment from the stack else return false
    cout<<"\nSolver started";
    env->var() = build_vars();

    size_t  var = 0;
    size_t new_line = 0;
    EqSolve<MX, V> solver;
    size_t nbrunches = 0;
    size_t nsolved = 0;
    //size_t brunch = 0;
    typename Env<MX, V >::ptr new_env;
    for(;;) {

        while(find_inequ(env->var(), env->line(), var)) { /// solve all inequlities
            new_line = env->line();
/*
            cout<<"\n<";
            cout<<"\ninq: "<<inq();
            cout<<"\nline, var: "<<new_line<<", "<<var;
            cout<<"\n>";
*/
            if (has_next(new_line, var)){
                new_env = typename Env<MX, V>::ptr(new Env<MX, V >(*env));
                new_env->line() = new_line;
                env_stack_.push(new_env);
//           cout<<"\nnew_line: "<<new_line;
            } else {
                /// nothing
            }
            
            solver.Solve(env->inq(), env->cnt(), env->line(), var);
            ++nsolved;
            
            
            env->var().erase(var);
            env->line() = 0;
        }
        /// return if proved

        if (check_solution())
            return InqSolverResult(true, nbrunches, nsolved);
        ++nbrunches;
        
//	cout<<"\rbrunch "<<nbrunches<<flush;
        if (!env_stack_.empty()) {
            env = env_stack_.top();
            env_stack_.pop();
        } else {
            break;
        }
    }
    /// no solutions
    return InqSolverResult(false, nbrunches, nsolved);
}


template <class MX, class V>
Vars InqSolverSubst<MX, V >::build_vars() const {
    Vars vars;
    /// Build set of varibles
    for(size_t j = 0; j<this->inq().size2(); ++j) {

#ifndef BOOST_MATRIX
        if (this->inq().begin1(j) != this->inq().end1(j)) 
#endif
#ifdef 	BOOST_MATRIX
	bool empty = true;
	typename MX::const_iterator1 it = (this->inq().begin2()+j).begin();
	typename MX::const_iterator1 it_end = (this->inq().begin2()+j).end();
	while( it != it_end && empty){
	    if (*it!= 0)
	    	empty = false;
	    	++it;
	}
	if (!empty) 
#endif
        {
            vars.insert(j);
        }
	

    }
   // cout<<"\nvariables builded: "<<vars;
    assert(vars.size()<=this->inq().size1());
    return vars;
}


/// Input:
///   varible (first in vars) -  which can be substituted by a several different ways
///   line - which is the starting position
/// Goal:
/// Find a next position to eliminate or return false if such position
/// does not exist.

template <class MX, class V>
bool InqSolverSubst<MX, V >::find_inequ( Vars & vars, size_t & line, size_t & var ) const {
//    cout<<"\nfind_inequ: \n"<<this->inq()<<flush;
    Vars::iterator it_var = vars.begin();
    typename MX::value_type zero = typename MX::value_type(0);
    while(it_var != vars.end()) {
        var = *it_var;
        bool found = false;
        while(line < this->inq().size1() && !found ) {
            if ((this->inq()(line, var)) != zero )
                found = true;
            else
                ++line;
        }

        if (found) {
            return true;
        }

        vars.erase(it_var);
        it_var = vars.begin();
        line = 0;
    }
    return false;
}

template <class MX, class V>
bool InqSolverSubst<MX, V>::check_solution( ) const {

#ifndef BOOST_MATRIX
    assert(this->inq().nlive()==0 && "Matrix shall be empty." );
#endif

#ifdef BOOST_MATRIX
    assert((nlive<MX>(this->inq())==0) && "Matrix shall be empty." );
#endif

    /// check solution
    typename V::const_iterator res = find_if(this->cnt().begin(), this->cnt().end(), std::bind2nd(std::greater<typename MX::value_type>(), typename MX::value_type(0)));
    if (res == this->cnt().end()) {
        return true;
    }
    return false;
}

/// ::has_next(size_t & line, const size_t & row)
/// For the given row, the function looking for the first nonzero element starting from
/// line+1 while end of the row is not reached
/// It does not metter how we choose a last varible, so this case was skiped

template <class MX, class V>
bool InqSolverSubst<MX, V >::has_next(size_t & line, const size_t & row ) const {

    if (row+1 < this->inq().size2()) { /// not a last varible

        typename MX::value_type zero(0);
        while(++line < this->inq().size1()) {
            if (this->inq()(line, row) != zero) {
                return true;
            }
        }
    }

    return false;
}

template <class MX, class V>
bool InqSolverSubstLEq<MX, V >::check_solution( ) const {
#ifndef BOOST_MATRIX
    assert(this->inq().nlive()==0 && "Matrix shall be empty." );
#endif

#ifdef BOOST_MATRIX
    assert((nlive<MX >(this->inq())==0) && "Matrix shall be empty." );
#endif

    /// check solution
    typename V::const_iterator res = find_if(
    	this->cnt().begin(), 
    	this->cnt().end(), 
    	std::bind2nd(
    		std::less<typename MX::value_type>(), 
    		typename MX::value_type(0)
    	)
    );
    
    if (res == this->cnt().end()) {
        return true;
    }
    return false;
}
